setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
options(stringsAsFactors = FALSE)
library(tidyverse)
library(cem)




### Due to privacy concerns, we cannot provide our raw data.
### This script replicates our matching procedure given
###    our list of cancer victims, opioid overdose victims, and their families.
### It then produces all the necessary plots.
### If you are interested in obtaining our raw data,
### Please email aaronkaufman@nyu.edu



###################################################################################################
###################################### Matching ################################################### 
###################################################################################################

# Load in the anonymized data with flags for opioid and cancer

load("cancer_families_anonymized.RData")

load("od_families_anonymized.RData")

load("ct_vf_13_17_anonymized.RData")
vot4$cancer = 0
vot4$cancer[vot4$id_voter %in% cancer$id_voter] = 1

vot4$opioid = 0
vot4$opioid[vot4$id_voter %in% od$id_voter] = 1

gdata::keep(vot4, od, cancer, sure=T)


# 2) Matching: cancer to full, opioid to full

vars = colnames(vot4)[colnames(vot4) %in% c("age", "sex", "imputed_race", "median_income.x",
                                            "county", "household_occupancy", "cd_party_13",
                                            "v_pres_general_12", "v_pres_primary_12",
                                            "v_cong_primary_10", "v_cong_general_10", "votprop", "cancer",
                                            "opioid", "id_voter")]

temp = vot4[,which(colnames(vot4) %in% vars)]
temp$cd_sex.x = as.factor(temp$cd_sex.x)
temp$v_pres_general_12 = as.factor(temp$v_pres_general_12)
temp$age[temp$age>100] = NA

m_od = cem(treatment = "opioid", data=temp[temp$cancer==0,-1])
m_cancer = cem(treatment = "cancer", data=temp[temp$opioid==0,-1])


## Did the matching work?
od_matched = vot4[m_od$matched,]
cancer_matched = vot4[m_cancer$matched,]
# Make sure that most of the observations are matched
sum(od_matched$opioid)/nrow(od)
sum(cancer_matched$cancer)/nrow(cancer)



###################################################################################################
###################################### Analysis ################################################### 
###################################################################################################

### Table 1

load("full_voterfile_anonymized.RData")
## Put them in a table, get p values

library(tidyverse)
library(xtable)

tab1 <- vot %>%
  group_by(indicator) %>%
  summarize(male = mean(cd_sex=="M"),
            white = mean(imputed_race=="white"),
            black = mean(imputed_race=="black"),
            hispanic = mean(imputed_race=="hispanic"),
            asian = mean(imputed_race=="asian"),
            age = mean(age, na.rm=T),
            dem13 = mean(cd_party_13=="D"),
            turnout_g12 = mean(v_pres_general_12))
tab1_overall  <- vot %>%
  summarize(male = mean(cd_sex=="M"),
            white = mean(imputed_race=="white"),
            black = mean(imputed_race=="black"),
            hispanic = mean(imputed_race=="hispanic"),
            asian = mean(imputed_race=="asian"),
            age = mean(age, na.rm=T),
            dem13 = mean(cd_party_13=="D"),
            turnout_g12 = mean(v_pres_general_12))

Ns = c(nrow(vot), sum(vot$indicator==1), sum(vot$indicator==2), sum(vot$indicator==3), sum(vot$indicator==4))

tab1[1,] = c(0,tab1_overall[1,])
tab1 = cbind(tab1, Ns)
colnames(tab1) = c("indicator", "Male", "White", "Black", "Hispanic", "Asian", "Age",
                   "\\% Democrat", "General Turnout 2012", "N")

tab1 = t(tab1)
colnames(tab1) = c("CT Pop", "Cancer Deaths", "Cancer Families", "OD Deaths", "OD Families")

writeLines(text=print(xtable(tab1[-1,], digits=3)), con="table1.tex")



# Plots
### 3a) Turnout

## I want to see avg change in partisanship from 2013 to 2017
## Among both the family and household dfs, and the overall voter file
## Same for turnout in general and primary 2012 and 2016

## hd is opioid, fd is cancer

# These two are the death families
cancer_a = cancer_matched[cancer_matched$cancer==1,]
opioid_a = od_matched[od_matched$opioid==1,]

# These two are the voterfile-matched samples
cancer_vf = cancer_matched[cancer_matched$cancer==0,]
opioid_vf = od_matched[od_matched$opioid==0,]

v = cancer_vf$v_pres_gen_16 == TRUE
w = opioid_vf$v_pres_gen_16 == TRUE 
x = cancer_a$v_pres_gen_16 == TRUE
y = opioid_a$v_pres_gen_16 == TRUE 
z = vot4$v_pres_gen_16 == TRUE 

vv = cancer_vf$v_pres_prim_16 == TRUE
ww = opioid_vf$v_pres_prim_16 == TRUE 
xx = cancer_a$v_pres_prim_16 == TRUE
yy = opioid_a$v_pres_prim_16 == TRUE 
zz = vot4$v_pres_prim_16 == TRUE 

#Pr(Vote16|R&opiod&Vote12) - Pr(Vote16|D&opiod&Vote12) > Pr(Vote16|R&cancer&Vote12)- Pr(vote16|D&cancer&Vote12)
a = cancer_a$v_pres_gen_16[cancer_a$cd_party_13=="D"] == TRUE # .725
b = cancer_a$v_pres_gen_16[cancer_a$cd_party_13=="R"] == TRUE # .827

c = opioid_a$v_pres_gen_16[opioid_a$cd_party_13=="D"] == TRUE # .747
d = opioid_a$v_pres_gen_16[opioid_a$cd_party_13=="R"] == TRUE # .803

e = vot4$v_pres_gen_16[vot4$cd_party_13=="D"] == TRUE # .804
f = vot4$v_pres_gen_16[vot4$cd_party_13=="R"] == TRUE # .830

g = cancer_vf$v_pres_gen_16[cancer_vf$cd_party_13=="D"] == TRUE # .725
h = cancer_vf$v_pres_gen_16[cancer_vf$cd_party_13=="R"] == TRUE # .827

i = opioid_vf$v_pres_gen_16[opioid_vf$cd_party_13=="D"] == TRUE # .747
k = opioid_vf$v_pres_gen_16[opioid_vf$cd_party_13=="R"] == TRUE # .803


aa = cancer_a$v_pres_prim_16[cancer_a$cd_party_13=="D"] == TRUE # .725
bb = cancer_a$v_pres_prim_16[cancer_a$cd_party_13=="R"] == TRUE # .827

cc = opioid_a$v_pres_prim_16[opioid_a$cd_party_13=="D"] == TRUE # .747
dd = opioid_a$v_pres_prim_16[opioid_a$cd_party_13=="R"] == TRUE # .803

ee = vot4$v_pres_prim_16[vot4$cd_party_13=="D"] == TRUE # .804
ff = vot4$v_pres_prim_16[vot4$cd_party_13=="R"] == TRUE # .830

gg = cancer_vf$v_pres_prim_16[cancer_vf$cd_party_13=="D"] == TRUE # .725
hh = cancer_vf$v_pres_prim_16[cancer_vf$cd_party_13=="R"] == TRUE # .827

ii = opioid_vf$v_pres_prim_16[opioid_vf$cd_party_13=="D"] == TRUE # .747
kk = opioid_vf$v_pres_prim_16[opioid_vf$cd_party_13=="R"] == TRUE # .803


#### Make Fig 1: Turnout
## three sets of points with error bars each
## Left: cancer D, voterfile-matched cancer D, cancer overall, voterfile-matched cancer overall, cancer R, voterfile-matched cancer R
## Middle: opioid D, voterfile-matched opioid D, opioid overall, voterfile-matched opioid overall, opioid R, voterfile-matched opioid R
## Right: Voterfile D, voterfile all, voterfile R
## color by party

points = sapply(list(a,g,x,v,h,b,  c,i,y,w,k,d, e,z,f), FUN=mean)
ses = sapply(list(a,g,x,v,h,b,  c,i,y,w,k,d, e,z,f), FUN=function(x){
  1.96*sd(x)/sqrt(length(x))
})

uppers = points + ses
lowers = points - ses

points2 = sapply(list(aa,gg,xx,vv,hh,bb,  cc,ii,yy,ww,kk,dd, ee,zz,ff), FUN=mean)
ses2 = sapply(list(aa,gg,xx,vv,hh,bb,  cc,ii,yy,ww,kk,dd, ee,zz,ff), FUN=function(x){
  1.96*sd(x)/sqrt(length(x))
})

uppers2 = points2 + ses2
lowers2 = points2 - ses2

pdf(file="paper/figures/fig1_v5.pdf", width = 7, height = 9)

par(mfrow=c(2,1), mai=c(0.4,.8,.4,.4))

plot(NA, xlim=c(0,3), ylim=c(.45,.75),ylab="2016 General Election Turnout", xaxt="n", xlab="")
points(x=c(.15,.2, .45, .5, .75,.8,  1.15, 1.2, 1.45, 1.5, 1.75,1.8, 2.15, 2.45, 2.75),
       y=points, pch=16, cex=1.5,
       col=c(rep(c("steelblue4", "steelblue1", "grey25", "grey70", "firebrick1", "firebrick4"),2),
             c("steelblue3", "grey", "firebrick3")))
segments(x0=c(.15,.2, .45, .5, .75,.8,  1.15, 1.2, 1.45, 1.5, 1.75,1.8, 2.15, 2.45, 2.75),
         x1 = c(.15,.2, .45, .5, .75,.8,  1.15, 1.2, 1.45, 1.5, 1.75,1.8, 2.15, 2.45, 2.75),
         y0 = lowers, y1 = uppers,
         col=c(rep(c("steelblue4", "steelblue1", "grey25", "grey70", "firebrick1", "firebrick4"),2),
               c("steelblue3", "grey", "firebrick3")), lwd=2)
abline(h=points[13], col="steelblue", lty="dashed", lwd=.5)
abline(h=points[15], col="tomato", lty="dashed", lwd=.5)
abline(h=points[14], col="grey", lty="dashed", lwd=.5)
abline(v=c(1,2), col="grey", lwd=0.5)
axis(side=1, at = c(.5, 1.5, 2.5), labels = c("Cancer", "Opioids", "Statewide"))
# legend
points(x = rep(2.55, 6), y=seq(.45,.525, length.out=6), pch=16, cex=1.5,
       col = rev(c("steelblue4", "steelblue1", "grey25", "grey70", "firebrick4", "firebrick1")))
text(x = rep(2.55, 6),
     y = seq(.45,.525, length.out=6), pos=4, cex=.6,
     labels = rev(c("Democrat Families", "Democrat Control",
                    "All Families", "All Control",
                    "Republican Families", "Republican Control")))
segments(x0 = 2.5, x1 = 2.5, y0 = 0, y1 = .54)
segments(x0 = 2.5, x1 = 3.5, y0 = .54, y1 = .54)

plot(NA, xlim=c(0,3), ylim=c(.1,.5),ylab="2016 Primary Election Turnout", xaxt="n", xlab="")
points(x=c(.15,.2, .45, .5, .75,.8,  1.15, 1.2, 1.45, 1.5, 1.75,1.8, 2.15, 2.45, 2.75),
       y=points2, pch=16, cex=1.5,
       col=c(rep(c("steelblue4", "steelblue1", "grey25", "grey70", "firebrick1", "firebrick4"),2),
             c("steelblue3", "grey", "firebrick3")))
segments(x0=c(.15,.2, .45, .5, .75,.8,  1.15, 1.2, 1.45, 1.5, 1.75,1.8, 2.15, 2.45, 2.75),
         x1 = c(.15,.2, .45, .5, .75,.8,  1.15, 1.2, 1.45, 1.5, 1.75,1.8, 2.15, 2.45, 2.75),
         y0 = lowers2, y1 = uppers2,
         col=c(rep(c("steelblue4", "steelblue1", "grey25", "grey70", "firebrick1", "firebrick4"),2),
               c("steelblue3", "grey", "firebrick3")), lwd=2)
abline(h=points2[13], col="steelblue", lty="dashed", lwd=.5)
abline(h=points2[15], col="tomato", lty="dashed", lwd=.5)
abline(h=points2[14], col="grey", lty="dashed", lwd=.5)
#segments(y0=points2[14], y1=points2[14], x0=-1, x1=2.55, col="grey", lty="dashed", lwd=0.5)
abline(v=c(1,2), col="grey", lwd=0.5)
axis(side=1, at = c(.5, 1.5, 2.5), labels = c("Cancer", "Opioids", "Statewide"))

points(x = rep(2.55, 6), y=seq(.1,.2, length.out=6), pch=16, cex=1.5,
       col = rev(c("steelblue4", "steelblue1", "grey25", "grey70", "firebrick4", "firebrick1")))
text(x = rep(2.55, 6),
     y = seq(.1,.2, length.out=6), pos=4, cex=.6,
     labels = rev(c("Democrat Families", "Democrat Control",
                    "All Families", "All Control",
                    "Republican Families", "Republican Control")))
segments(x0 = 2.5, x1 = 2.5, y0 = 0, y1 = .22)
segments(x0 = 2.5, x1 = 3.5, y0 = .22, y1 = .22)


dev.off()







### 4) Party change


# Two graphs: turnout (include other elections too), PID (transition matrix? fancier thing?)
## The point of this graph, like Fig 1, is the compare relative changes
## The point is that Dems defect basically the same amount, but Rs both defect more
## So I'd need 3 transition networks side by side? Seems complicated


## Compress Party ID to R, D, or I
vot4$cd_party_13[!vot4$cd_party_13 %in% c("R", "D")] = "I"
vot4$cd_party_17[!vot4$cd_party_17 %in% c("R", "D")] = "I"
opioid_a$cd_party_13[!opioid_a$cd_party_13 %in% c("R", "D")] = "I"
opioid_a$cd_party_17[!opioid_a$cd_party_17 %in% c("R", "D")] = "I"
cancer_a$cd_party_13[!cancer_a$cd_party_13 %in% c("R", "D")] = "I"
cancer_a$cd_party_17[!cancer_a$cd_party_17 %in% c("R", "D")] = "I"
opioid_vf$cd_party_13[!opioid_vf$cd_party_13 %in% c("R", "D")] = "I"
opioid_vf$cd_party_17[!opioid_vf$cd_party_17 %in% c("R", "D")] = "I"
cancer_vf$cd_party_13[!cancer_vf$cd_party_13 %in% c("R", "D")] = "I"
cancer_vf$cd_party_17[!cancer_vf$cd_party_17 %in% c("R", "D")] = "I"




fullmat = matrix(nrow=3, ncol=3, byrow=F,
                 c(mean(vot4$cd_party_17[vot4$cd_party_13 == "R"] == "R", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "R"] == "I", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "R"] == "D", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "I"] == "R", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "I"] == "I", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "I"] == "D", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "D"] == "R", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "D"] == "I", na.rm=T),
                   mean(vot4$cd_party_17[vot4$cd_party_13 == "D"] == "D", na.rm=T)))
colnames(fullmat) = rownames(fullmat) = c("R", "I", "D")
fullmat = round(fullmat, 3)


cancermat = matrix(nrow=3, ncol=3, byrow=F,
                   c(mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "R"] == "R", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "R"] == "I", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "R"] == "D", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "I"] == "R", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "I"] == "I", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "I"] == "D", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "D"] == "R", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "D"] == "I", na.rm=T),
                     mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "D"] == "D", na.rm=T)))
colnames(cancermat) = rownames(cancermat) = c("R", "I", "D")
cancermat = round(cancermat, 3)

odmat = matrix(nrow=3, ncol=3, byrow=F,
               c(mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "R"] == "R", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "R"] == "I", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "R"] == "D", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "I"] == "R", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "I"] == "I", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "I"] == "D", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "D"] == "R", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "D"] == "I", na.rm=T),
                 mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "D"] == "D", na.rm=T)))
colnames(odmat) = rownames(odmat) = c("R", "I", "D")
odmat = round(odmat, 3)

cancermat_vf = matrix(nrow=3, ncol=3, byrow=F,
                   c(mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "R"] == "R", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "R"] == "I", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "R"] == "D", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "I"] == "R", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "I"] == "I", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "I"] == "D", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "D"] == "R", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "D"] == "I", na.rm=T),
                     mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "D"] == "D", na.rm=T)))
colnames(cancermat_vf) = rownames(cancermat_vf) = c("R", "I", "D")
cancermat_vf = round(cancermat_vf, 3)

odmat_vf = matrix(nrow=3, ncol=3, byrow=F,
               c(mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "R"] == "R", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "R"] == "I", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "R"] == "D", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "I"] == "R", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "I"] == "I", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "I"] == "D", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "D"] == "R", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "D"] == "I", na.rm=T),
                 mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "D"] == "D", na.rm=T)))
colnames(odmat_vf) = rownames(odmat_vf) = c("R", "I", "D")
odmat_vf = round(odmat_vf, 3)







tab2 = cbind(c(fullmat), c(cancermat), c(cancermat_vf), c(odmat), c(odmat_vf))
colnames(tab2) = c("All CT", "Cancer Families", "Cancer Matched", "OD Families", "OD Matched")

r_notr = c(mean(vot4$cd_party_17[vot4$cd_party_13 == "R"] != "R", na.rm=T),
           mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "R"] != "R", na.rm=T),
           mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "R"] != "R", na.rm=T),
           mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "R"] != "R", na.rm=T),
           mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "R"] != "R", na.rm=T))
r_notr = round(r_notr, 3)


d_notd = c(mean(vot4$cd_party_17[vot4$cd_party_13 == "D"] != "D", na.rm=T),
           mean(cancer_a$cd_party_17[cancer_a$cd_party_13 == "D"] != "D", na.rm=T),
           mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 == "D"] != "D", na.rm=T),
           mean(opioid_a$cd_party_17[opioid_a$cd_party_13 == "D"] != "D", na.rm=T),
           mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 == "D"] != "D", na.rm=T))
d_notd = round(d_notd, 3)

notr_r = c(mean(vot4$cd_party_17[vot4$cd_party_13 != "R"] == "R", na.rm=T),
           mean(cancer_a$cd_party_17[cancer_a$cd_party_13 != "R"] == "R", na.rm=T),
           mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 != "R"] == "R", na.rm=T),
           mean(opioid_a$cd_party_17[opioid_a$cd_party_13 != "R"] == "R", na.rm=T),
           mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 != "R"] == "R", na.rm=T))
notr_r = round(notr_r, 3)


notd_d = c(mean(vot4$cd_party_17[vot4$cd_party_13 != "D"] == "D", na.rm=T),
           mean(cancer_a$cd_party_17[cancer_a$cd_party_13 != "D"] == "D", na.rm=T),
           mean(cancer_vf$cd_party_17[cancer_vf$cd_party_13 != "D"]== "D", na.rm=T),
           mean(opioid_a$cd_party_17[opioid_a$cd_party_13 != "D"] == "D", na.rm=T),
           mean(opioid_vf$cd_party_17[opioid_vf$cd_party_13 != "D"] == "D", na.rm=T))
notd_d = round(notd_d, 3)

tab2 = rbind(r_notr, d_notd, notr_r, notd_d, tab2)
rownames(tab2) = c("R asd !R", "D asd !D", "!R asd R", "!D asd D", "R asd R", "R asd I", "R asd D",
                   "I asd R", "I asd I", "I asd D", "D asd R", "D asd I", "D asd D")

## Add tx and diff-in-diff
tab2 = as.data.frame(tab2)
tab2$cancer_tx = tab2[,2] - tab2[,3]
tab2$od_tx = tab2[,4] - tab2[,5]
tab2$difindif = tab2$od_tx - tab2$cancer_tx
tab2$gap = ""

tab2 = tab2[,c(1,2,3,6,9,4,5,7,8)]

tab2 = print(xtable::xtable(tab2, digits=3)) # these two lines are a cheap hack to get xtable to output $\to$
tab2 = gsub("asd", "$\\\\to$", tab2)

writeLines(text=tab2, con="./paper/tables/table2_revision.tex")


## Hypothesis tests
vot4$party_code = paste0(vot4$cd_party_13, vot4$cd_party_17)
opioid_a$party_code = paste0(opioid_a$cd_party_13, opioid_a$cd_party_17)
cancer_a$party_code = paste0(cancer_a$cd_party_13, cancer_a$cd_party_17)
opioid_vf$party_code = paste0(opioid_vf$cd_party_13, opioid_vf$cd_party_17)
cancer_vf$party_code = paste0(cancer_vf$cd_party_13, cancer_vf$cd_party_17)

vot4$category = "statewide"
opioid_a$category = "opioid"
cancer_a$category = "cancer"
opioid_vf$category = "opioid_matched"
cancer_vf$category = "cancer_matched"

vot5 = rbind(vot4[,c("party_code", "category")],
             opioid_a[,c("party_code", "category")],
             cancer_a[,c("party_code", "category")],
             opioid_vf[,c("party_code", "category")],
             cancer_vf[,c("party_code", "category")])

chisq.test(table(vot5$party_code[vot5$category %in% c("opioid", "cancer")],
                 vot5$category[vot5$category %in% c("opioid", "cancer")]))

chisq.test(table(vot5$party_code[vot5$category %in% c("opioid", "statewide")],
                 vot5$category[vot5$category %in% c("opioid", "statewide")]))

chisq.test(table(vot5$party_code[vot5$category %in% c("statewide", "cancer")],
                 vot5$category[vot5$category %in% c("statewide", "cancer")]))

chisq.test(table(vot5$party_code[vot5$category %in% c("cancer", "cancer_matched")],
                 vot5$category[vot5$category %in% c("cancer", "cancer_matched")]))

chisq.test(table(vot5$party_code[vot5$category %in% c("opioid", "opioid_matched")],
                 vot5$category[vot5$category %in% c("opioid", "opioid_matched")]))







library(GGally)
library(network)
library(sna)
library(ggplot2)
library(igraph)

bip = data.frame(D = c(0,1,1),
                 I = c(1,0,1),
                 R = c(1,1,0))

net=graph.adjacency(as.matrix(bip),mode="directed",weighted=TRUE,diag=T) 
V(net)$color = c("steelblue", "grey", "tomato")
E(net)$label = c(paste("Full: ", fullmat[2,1], "\nCancer: ", cancermat[2,1], "\nOpioids: ", odmat[2,1], sep=""),
                 paste("Full: ", fullmat[3,1], "\nCancer: ", cancermat[3,1], "\nOpioids: ", odmat[3,1], sep=""),
                 paste("Full: ", fullmat[1,2], "\nCancer: ", cancermat[1,2], "\nOpioids: ", odmat[1,2], sep=""),
                 paste("Full: ", fullmat[3,2], "\nCancer: ", cancermat[3,2], "\nOpioids: ", odmat[3,2], sep=""),
                 paste("Full: ", fullmat[1,3], "\nCancer: ", cancermat[1,3], "\nOpioids: ", odmat[1,3], sep=""),
                 paste("Full: ", fullmat[2,3], "\nCancer: ", cancermat[2,3], "\nOpioids: ", odmat[2,3], sep=""))

pdf(file="./paper/figures/fig2.pdf", height=9, width=15/2)
plot.igraph(net, edge.curved=TRUE, arrow.size=.3,
            layout=layout.fruchterman.reingold, vertex.size=24,
            edge.label.cex=1.5, vertex.label.cex=1.5)
dev.off()

